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1. Introduction 

Contemporary simulations of full QCD on the lattice typically use two degenerate light 
flavors ('up' and 'down') and one heavier flavor ('strange') of sea quarks. For each flavor the 
fermionic determinant is replaced by an integral over bosonic fields (f), the pseudofermions [1] 



detLioc / [#t][d0]e-<^^i<^ . (1.1) 



This identity requires all eigenvalues of the matrix D to have a positive real part. 

Most lattice Dirac operators obey 75-Hermiticity, D'^ = 75-D75, and so detD is real. 
However, the eigenvalues of lattice Dirac operators are typically complex and their real 
parts may not be positive-definite. Then the exponential in Eq. (1.1) cannot be interpreted 



as a conventional probability measure. This is the case for Wilson fermion actions. They 
have eigenvalues that are paired complex-conjugates or unpaired and real: 



detL»= JJ(|Ajf + m2)JJ(Ar + m) . (1.2) 



pairs 



The absence of chiral symmetry means that the sign of the real eigenvalues is not protected, 
so the determinant can have either sign. The solution to this problem is to simulate two 
degenerate fiavors at a time, that is, to rewrite their fermion determinant det D = det D^D 
where D^^ D fulfills the requirements of Eq. 1.1. So Eq. 11.21 becomes 



deiD^D = J[{\\i\^ + m^)W{\r+mf . (1.3) 

i r 

where the product index i runs over all modes. Simulations of a single flavor are, in general, 
not as straightforward, and one is forced to use either the Refreshed Molecular Dynamics 
or R algorithm [2], which is not exact, or the Rational Hybrid Monte Carlo (RHMC) 



algorithm [3], which uses the square root of the two flavor operator, Q = vW^D. This 
could introduce a systematic error if Xr + m could change sign. 

For the overlap Dirac operator [4, 5] there is a more natural solution. It exploits the 
chiral properties of the formulation, i.e. that the eigenmodes of D^D can be chosen to be 
chiral, with an eigenmode of each chirality per eigenvalue |Ajp + m^. A single flavor of 
overlap fermions can be simulated with one chiral pseudofermion. One only has to correct 
for the effect of the real modes. The real spectrum is known if one knows the topological 
charge of the configuration as defined by the index of the Dirac operator. The idea is based 
on a suggestion which can already be found in Ref. [6] but shall now be formulated and 
investigated more precisely. 

A simulation of any number rij of flavors thus involves a set of Uf chiral pseudo- 
fermion fields. It is also necessary to keep track of the global topology of the configuration 
during the simulation. Besides allowing simulations for any Uf, the chiral algorithm has an 
advantage over the traditional "non-chiral" Hybrid Monte Carlo (HMC) algorithm (with 
a Dirac spinor pseudo-fermion for each degenerate flavor pair): Away from topological 
boundaries, or in simulations in which topological changes are forbidden, it can eliminate 
part of the critical slowing down in the chiral limit and stabilize the inversion in sectors of 
non-trivial topology. This is achieved simply by running in the sector of chirality which is 
topologically trivial. 

In Sec. 2 we describe the method. The new ingredient, not mentioned in Ref. [6], is 
a technique for generating the initial chiral pseudofermion. In Sec. 3 we present results of 
some test simulations in small volumes. Some of these results are surprising, so in Sec. 4 
we describe a little solvable model which reproduces these features. 

2. The method 

Let us first fix our conventions and repeat a few properties of Neuberger's overlap operator. 
It is given by 

D = Do,{m = 0) = i?o [1 + l^e{h{-Ro))\ (2.1) 

with e{h) = h/vh? the sign function of the Hermitian kernel operator h = 75^ which is 
taken at the negative mass Rq. Its spectrum is symmetric, i.e. each non-real eigenvalue A is 
paired with its complex conjugate A*. The modes at zero and 2Rq are the only real modes. 
They are chiral and the excess of the zero modes of negative chirality over the ones with 
positive chirality gives the topological charge Q. Experience shows that there are always 
only zero modes of one chirality. For simplicity of the argument we will assume this in the 
following even though this assumption is not necessary. The squared Hermitian overlap 
operator H^ = (75!))^ = D'^ D commutes with 75 and therefore can have eigenvectors with 
definite chirality. The modes at zero and 4i?Q aside, the spectrum is doubled with a positive 
and a negative chirality eigenvector to each eigenvalue |Ap. It is therefore convenient to 
define the chiral projections {P± = 2(1^75)) so that the massive squared Hermitian overlap 
operator, with the usual convention for the mass terms, is 

2 
Hl{m) = P±H\m)P± = 2{rI - ^)P±(1 + e{h))P± + n?P± . (2.2) 



Let us call the chiral sector without the zero modes the sector of opposite chirality and 

opp- 



the associated Hermitian Dirac operator H^ . The chiral sector with the chirality of the 



zero modes shall be called the same chirality sector and -f^^ame the associated operator. 

The method which we are going to discuss exploits the fact that H^ has the same 
spectrum in both chiral sectors. Only the modes at m? and 4i?Q differ. The fermion 
determinant for one flavor is therefore the determinant of H^ in one chiral sector times a 
correction factor for the modes at m and 2Rq. This is summarized in the following relation 

detD = (m/2i?o)l^ldeti:f2^p = {2R^/m)\<^\deiHl^^ . (2.3) 

Since -ffopp ^'^d -f^^amc 9-^6 positive operators, it is thus straightforward to simulate a single 
flavor in hybrid Monte Carlo. 

First, however, one has to decide which action to simulate, i.e. how to choose the 
optimal chirality for the pseudofermion a given gauge configuration. To begin, for any 
configuration with topological charge, the inversion of the Dirac operator is computationally 
cheaper in the sector of opposite chirality. Typically, the non-zero modes are repelled by the 
zero mode(s) and the smallest non-zero mode is significantly above zero. The conditioning 
number of the Dirac operator (the ratio of largest to smallest eigenvalue) is therefore 
lower in the opposite chirality sector because the largest eigenvalue is always close to 2Rq 
regardless of chirality. (For H{rn)^ it is proportional to l/{m?' + A^j„).) 

During the molecular dynamics evolution the topological charge can change. One has 
to decide what to do when the chirality of the zero-modes flips (e.g. starting in a trajectory 
in a configuration with Q = \ and evolving through a region with Q = Q into Q = —1). 
There are in principle two options: One can keep the chirality of the operator to be fixed, 
or one can choose the chirality of the operator to be in the opposite chirality sector of the 
configuration. This means that when one changes the topology of the configuration, one 
also changes the chirality of the pseudofermion. The associated action is given by 

det D = {m/2Ro)\Q\detH^^p cc j m[d4>^] exp(|Q| log ^ - <^^i7opp0) • (2-4) 

To simulate this action we use the Hybrid Monte Carlo algorithm. We build on experi- 
ences previously published in Refs. [7-9]. The initial formulation of the HMC algorithm [10] 
for the overlap operator is given in Ref. [11]. At the start of a trajectory, one has to perform 
a heat-bath for the pseudo-fermion fields. In a conventional n/ = 2 simulation, one would 
cast a random vector ^ in both chiralities and compute (j) = H^. Then one would break (j) 
into its separate chiralities and use 

detD^ = /'[(i(/)+][d(/.^][#_][#l]exp(-(/)^i7+2^+ - (/.Ifi^-)- (2.5) 



However, we need a set of chiral (j) fields, of chirality o", chosen by heat bath. (How to 
do this is not described explicitly in Ref. [6].) We achieve this goal by generating a set of 
chiral Gaussian random fields ^a- We then define (pa = ^Ju^^a- To construct (\>cf we use a 
rational approximation to the square root in the region [m^,4i?Q], i.e. we approximate it 



by 

I ^'^^^i 

where we use the 6; and c; from the Zolotarov approximation to the sign function (there 
used to compute e{x) = x/Vx^). 

When we use Hasenbusch preconditioning [12], the action for the lower mass pseud- 
ofermions is 

Sf = ^IffT^^P^- (2-7) 



Our chiral pseudofermion is then taken to be y^ Ha{m)'^ / Hfj{m')'^£,ij, which we again ap- 
proximate by a Zolotarov formula, whose independent variable x obeys 

x-^ = H^{m'f/H^{mf = a + f3/H^{mf. (2.8) 

The range of x is {m/m')'^ to 1. 

During the trajectory one has to deal with the discontinuity due to the sign function 
in the definition of the overlap operator, when an eigenmode |Ao) of the kernel operator 
h{—Ro) changes sign. Fodor et al. [11] proposed a method of how to deal with this problem. 
One measures the height of the step in the effective action and then reflects or refracts the 
gauge field momentum as in classical mechanics. The computation of the height of the step 
is a major part of the total cost of the simulation. If we use the same chirality of Hf^ on 
both sides of that boundary this amounts to the change 

H^im) -^ Hl{m) ± {ARI - m^)P„\\^){\o\P^ = H^{m) . (2.9) 

From the Sherman-Morrison formula, 

1 1 /S^ 1 1 

Pa\Xo){\o\Pa^pTr^, (2.10) 



so the height of the step is given by 



A 



'^'^^.(m)^^'^' 



Finally, for completeness, the exact ratio of determinants is 

deiHlim] 



^^ ""^-TT^^-I^o)!'. (2.11) 



deiHlim) 



1 + 6CL. (2.12) 



In the abbreviated formulas 6 is the sign in Eq. p.9| , C = (4i?Q — m"^) and L is the matrix 
element {Xo\PtTH~'^{m)Pa-\Xo). 

Eq. |2.1l| is obviously only applicable if we use the same chiral sector on both sides. 
Otherwise, one has to run the inversion twice, which is very expensive and numerically less 
under control. This will occur during crossings into Q = 0. Technically, a change in the 
chirality which we are using for H^ amounts to a change in the chirality which we use for 



the embedding of the two-component chiral source into the four-component Wilson vector. 
We will therefore use the term source chirality for chirality of H^ . 

Since choosing the chirality which we use in the Q = Q sector depending on the 
configuration we start the trajectory from would violate reversibility, we choose it randomly 
at the beginning of each trajectory. So at a topological boundary, where we propose a 
tunneling into (5 = 0, the new source chirality will be different from the initial chirality 
half the time. The AS* which results from flipping the source chirality is very large and 
most of these crossings will be rejected. This suggests two other algorithms for simulating 
any number of flavors: 

• Simply restrict the simulation to a particular topological sector. There are simu- 
lational situations where this restriction is desirable. They include the calculation 
of the condensate using random matrix theory, or calculations of full QCD in the 
so-called epsilon regime. It is unknown whether these simulations are ergodic. If the 
manifold of gauge fields corresponding to fixed topology were smoothly connected, 
then HMC would (in principle) carry us from any gauge configuration to any other 
one by a series of small steps. However, if sectors of fixed topology were disjoint, or 
could be connected only by passage a sector of some other topological charge, HMC 
in a sector of fixed topology would not be ergodic. We are aware of no proofs one 
way or the other, for four dimensions. Arguments we construct based on instanton 
phenomenology, where Q counts the excess of instantons over anti-instantons, argue 
that there is no problem: in the different configurations, the location of the odd in- 
stanton(s) moves around, and their sizes shrink and grow, but this is all continuous. 
So is the appearance of pairs of instantons and anti-instantons, as they grow from 
fluctuations of a single plaquette, or annihilate similarly. To produce the entire func- 
tional integral, results from different topological sectors can be combined using Eq. 
p. 12 , as described by Ref. [13]. 



Alter the tunneling probability and reweight the resulting data set, if necessary. A 
simple way to do this is to pick all chiral sources to carry the same chirality and 
begin the simulation either in Q = or in a topological sector in which the sources 
do not have zero modes. Here we must assume that configurations with zero modes 
in both chiralities never appear. Then allow topological changes in which the sources 
do not have zero modes, but prohibit transitions which would create zero modes in 
the source chirality. For example, we could set the source chiralities to be positive 
and only allow transitions into Q = n_ — n+ > 0. (We will call this the "fixed 
chirality algorithm.") Unless there are disconnected sectors at Q 7^ which can 
only be reached by some passage through Q < (for example Q = l— >Q = 0— > 
(5 = — 1— >(5 = 0— >Q = 1— >(5 = 2) the algorithm will generate an ensemble 
with the correct Boltzmann weighting ratio between sectors of all Q > Q. Under a 
parity transformation a gauge configuration with positive Q is converted into one 
with negative Q. In the analysis of an ensemble generated with this algorithm, 
measurements on the Q = ^ configurations need to be reweighted with a factor 1/2 
compared to those from configurations with non-trivial topology. 



At a topological boundary we add A\Q\ log(m/(2i?o)) to the pseudofermion AS* before 
deciding whether to reflect or refract. As a variant on this approach (which we have 
not tried) one could leave out a fraction of the |(5| log(m/(2i?o)) factor from the 
action during the HMC evolution and include it later with a real reweighting. 

3. QCD simulations 

We make tests on four different data sets: Set A is generated with the standard two-flavor 
HMC algorithm on an 8'^ x 6 lattice, at a lattice spacing of a ~ 0. 16 fm and a bare quark mass 
of am = 0.05. This is close to the crossover between the chirally broken low temperature 
phase and the chirally restored high temperature phase. The details of this simulation 
and its parameters are similar to those discussed in Ref. [8]. We only mention that our 
gauge connections are stout links [14] and that we use Hasenbusch preconditioning [12] 
with one extra pair of pseudo-fermion fields at a higher mass. Here we use three levels of 
stout smearing. The molecular dynamics integration has a multiple time step integration 
similar to that of Ref. [15]. Set A' is run at the same parameters but using the new chiral 
algorithm. Set B is generated with the fixed chirality algorithm. We restrict these studies 
to Uf = 2 because we can make comparisons to the usual (nonchiral) HMC algorithm. (We 
have also done extensive running with n/ = 1, which we will report elsewhere.) 

Finally, we have run the new algorithm to generate sets (labeled C) of 10*^ lattices at the 
same lattice spacing as sets A and A', with two steps of stout smearing, at fermion masses 
of amq = 0.05, 0.03 and 0.015. These simulations are done in sectors of fixed topology 
by switching off the possibility of refraction. This allows us to study the behavior in the 
different topological sectors. Otherwise, in particular at small quark masses, the fermion 
determinant suppresses the sectors of non-zero topology and it is hard to get sufficient 
statistics there. 

We check that our starting pseudofermion field is chosen appropriately by computing 
(f)aHa{m)^'^(l)a- (or whcu we use Hasenbusch preconditioning, (l)]yH'^{rn')Ha{rn)^'^(j)a-) and 
comparing this value to the heat bath initialization ^Xia- With high accuracy evaluations 
of the Zolotarov formula and a tight convergence criterion for the Conjugate Gradient 
inversion of the quark propagator [r^r = 10~^^), the deviation in the action from its heat 
bath value is held below 0.02 or so, out of a total fermion energy in our simulations of a 
few times 10^. This is small compared to the typical violation of energy conservation in 
our molecular dynamics trajectory. 

Let us turn to the critical slowing down and the cost of the inversion (we do not 
make statements about the auto-correlation time because our data set is too small to make 
a definite statement). The new method has an overhead at the start of the trajectory 
because one has to hit the Gaussian source with the square root, which involves a multi- 
mass inversion of the overlap operator, instead of just the operator H. However, this is 
only a small fraction of the total cost of the algorithm. 

Since the simulation for set C is done in a fixed topological sector, we have significant 
statistics for the Q = ±1 sector for smaller quark masses. We can thus study the critical 
slowing down of the inversion at trivial and non-trivial topology. For all three quark masses 



we found no significant difference in the cost of the inversion between the Q = and the 
Q = ±1 runs. At fixed topology, the variation of the number of CG steps for the three 
different bare quark masses {am = 0.015, 0.03 and 0.05) is below 20%. We can conclude 
(at least for the small volumes and for our parameters) that critical slowing down is largely 
eliminated by the chiral algorithm. 

The other algorithms allow for topological changes. We made simulation runs of 250- 
300 trajectories for each algorithm. We observed that the plaquettes from all simulations 
were consistent within uncertainties. In all three runs we had tunnels from Q = into and 
out of IQI = 1. The tunneling rate was too low in all three simulations to say anything 
meaningful about autocorrelation times. All three runs used the same parameters. All had 
acceptance rates of about 90 per cent. All had about 1.8 attempted topological changes 
per trajectory. The cost of a trajectory in units of the number of applications of H^ to a 
trial vector were 2430(23) for set A, 2798(47) for A', and 2221(17) for B. The excess of 
A' over A is due to the startup. The decrease of B from A is from the modest decrease in 
the conditioning number because we never run in a topologically nontrivial sector. 

The algorithm refracts when AS, the change in action, is smaller than ^{N\tt)'^, the 
squared projection of the gauge momentum normal to the surface of topology change. We 
first show a histogram of (A^JTr)^ in Fig. |^. As expected there is little difference in this 
quantity between the algorithms. 

Next we look at AS and show histograms of this quantity in Fig. 0. These distributions 
are quite different. To make sense of them, we realize that while for the nonchiral algorithm 
all crossings are "similar," in the sense that all transitions involve contributions from zero 
modes in either the initial or final state, that is not the case for the chiral algorithms: 
either the transition involves a change in topology in which the zero mode has appeared or 
disappeared from the opposite chirality sector, or in the same sector as the simulation. In 
the case of an opposite chirality change, the magnitude of the topology can either increase 
or decrease. 

In Fig. m we further divide the contributions to Fig. ^d and show histograms of AS" 
from algorithm A' for two cases. In panel (a) we show AS from transitions when the 
running chirality is different from the topology of the proposed crossing. This histogram 
itself has two components, a narrow one and a wide one, which we will shortly separate. 

In panel (b) we show AS* for transitions in which the final topology and the running 
topology have the same sign. In this algorithm we compute the action in the new sector 
by flipping the pseudofermion chirality. This gives a very noisy estimator for AS, with a 
large mean and deviation. 

Changes "up" (Q = ^ Q = 1) and "down" (Q = 1 ^ Q = 0) for the fixed 
chirality algorithm (data set B) are illustrated in Fig. ^ These distributions are also quite 
asymmetric. A breakdown of Fig. ^ would duplicate this figure. 

The common feature of these plots is that the distribution of AS is large and wide 
when the lowest eigenvalue of H^ would shrink if the topology changed, and is small and 
narrow when the eigenvalue would grow. The lowest eigenvalue in the same-chirality sector 
shrinks when a zero mode appears. In the opposite-chirality sector, when the magnitude 
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Figure 1: Histograms of (A^JTr)^, the squared normal component of the gauge momentum at a 
topological boundary, from the non-chiral algorithm A (a), the chiral algorithm A' (b), and the 
fixed chirality algorithm B (c). 



of the topological charge increases the smallest eigenvalue also increases, and when the 
magnitude of |(5| drops, so does the smallest eigenvalue. 

This is most apparent in Fig. |^. The distribution is wide in Fig. ^ because when 
we attempt to tunnel out of Q = 0, a near-zero mode appears in the spectrum. Fig. ^a 
contains many low AS values (Q = to Q = 1) and a few high values (Q = 1 to Q = 0) 
as in Fig. |^. The two components in Fig. gb are a narrow one for transitions from Q = 
to Q = 1 with Q = —1 sources and a wider one, one for Q = 1 to Q = transitions. This 
second component has a contribution in which the Q = sector's chirality is flipped. We 
don't have enough data from the chiral algorithms for transitions from Q = 1 to Q = 2 to 
make a histogram, but there is a strong hint of a small AS" for transitions from 1 to 2 and 
a big AS for transitions down. 

In all of these distributions, the size of the fluctuations in AS for a particular kind of 
proposed topological change is strongly correlated to the size of AS" itself. This correlation 
arises because AS comes from an average over a set of Gaussian random vectors and 
because the change in the pseudofermion action is limited to a single crossing mode. (We 
are also assuming that the initial and final pseudofermion chiralities are identical.) Taking 
one pseudofermion and assuming that we have refreshed immediately before encountering 
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Figure 2: Histograms of AS, the the height of the step at a topological boundary, from the 
non-chiral data set A (a), the chiral data set A' (b), and the fixed chirality data set B (c). 
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Figure 3: Histograms of AS from algorithm A', where the proposed topology change is from (a) 
transitions away from the running chirality and (b) into the running chirality. 



a crossing, the action change per flavor is 



1 



AS = C{^/H!g^VH!-'^)C^ 



(3.1) 
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Figure 4: Histograms of AS" from data set B, where the proposed topology change is from (a) 
to 1 and (b) 1 to 0. 



where l/-ff^ is given by Eq. ^^. The average is done with respect to a weight factor which 
is a pure Gaussian, 

fd^U( O exp(-e^O 



(O) 



Jd^d^eM-^^O ' 



(3.2) 



so if O ~ C V^j then (O) = TvV. The trace runs over one state, the crossing state, and so 
TVy is the number Vq or 

Because the measure is just a Gaussian, the squared variance is o"^ = ((AS*)^) — (A5)^ = 
Try^. Again, only the crossing state contributes to V, TrV^ = Vq, and so the variance a 
is equal to the absolute value of ( A5) . 

Because a topological boundary can only be crossed when {N-H)'^ > 2 AS, and because 
{N ■H)'^ is always on the order of unity (recall Fig. ||), only when AS is in its low- value tail 
can a crossing occur. However, because the average value of AS* is equal to its fluctuation, 
this is a constant fraction of the AS* sample. This is how the algorithm preserves detailed 
balance. 



4. A model calculation 

To illustrate our results, we have constructed a solvable model with a discontinuity in its 
spectrum. It is a simple system which is confined to a box and inside of that box has two 
regions with different weights. To be specific, the partition function is 



dx < 



oo 



oo 



X < -1 



detM| -1 < X < 
det M| < X < 1 



X > 1 



(4.1) 



10 



For X < 0, the weight is given by the determinant of 




and for x > simply by 



Mr 




(4.2) 



(4.3) 



Besides having different determinants, the matrices do not commute, so their eigenvectors 
change across the step at x = 0. Since Ml has eigenvalues m and 2 + tti it plays the role 
of the sector of QCD with the lower eigenvalue of the Dirac operator. 

We simulate this theory with the HMC algorithm of Ref. [11], reflecting off the walls 
and reflecting/refracting at the step at x = 0. We introduce the determinant (s) by pseudo- 
fermions 

(detMi)2= /"#d(/.'^exp(-(/.t(M/Mi)~V)- (4.4) 

We allow for the possibility of Hasenbusch preconditioning. For n pseudofermion masses, 
the jth mass is rrij = m^"-~^+^ii^. 

The "gauge field" is the position variable x which we drive with a momentum p. We 
want to see two things: (i) What is the width of AS" in barrier crossings? (ii) What affects 
the tunneling rate? Fig. |5| shows AS at the crossing from a simulation in the model 
with ?TT, = 0.1 and a single pseudofermion. The agreement with what we saw in the QCD 
simulation - a wide distribution when the minimum eigenvalue drops, a narrow distribution 
when it rises - is striking. 
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Figure 5: AS at the crossing in the model, with m = 0.1 and a single pseudofermion. 

Now for the tunneling rate: we first compute the refraction probability in this model. 
In the exact case, the rate of tunnels from right to left is 



/"OO 

P{R -^L)= Nlr / pdpexp{-p^)9ip^ - \ogidet{Mn/MLf)) 
Jo 



(4.5) 



11 



where Nm is the number of left moving particles on the right side and the extra factor 
of p in the integrand counts the flux across the barrier. The tunneling rate in the other 
direction is 

r-oo 

P{L -^R) = Nrl / pdpeM-P^) (4.6) 

where Nm is the number of right movers on the left side. The theta-function is absent 
because the logarithm is always negative, and the integral is just unity: the crossing rate 
is 100 per cent. Of course, Nlr = Nrb, and Nrl = ^ll from reflections at the ends. 
Equating the tunneling rates gives Nm = A'^rx det(Mi/M/j)^ which is the statement of 
detailed balance. 



Eq. 4.6 represents an upper bound on the tunneling rate from right to left, and the 
tunneling rate for HMC will always be less than this value. Since detailed balance is obeyed, 
the tunneling rate in the other direction is also suppressed. 

In the stochastic case and if we are on the left side, we choose (f) = M^^ where ^ is 
Gaussian. Then AS" = i'^{W - 1)^ where W = Ml{m\^Mr)-^mI. We can rotate the 
pseudofermion integration measure to a basis which diagonalizes W and the tunneling rate 
is 

poo p' 

P{L ^R) = Nrl / pdpexp( V) H ^^^ «^P(- E ^*')^(p' " E ^^^^^ " l))' (4.7) 

If the eigenvalues of W were all less than unity, the momentum integral would be uncon- 
strained. That does not happen, however: it is easy to show that the eigenvalues oiW — 1 
are ei = (3 + 2m)/(l + m)^ and e2 = —(1 + 2m)/(l + rn)^ . The tunneling rate is reduced 
from Nlr to Nlr[1 — ef /(I + ei)(ei — €2)]- To preserve detailed balance, the tunneling rate 
in the opposite direction must be suppressed by the same amount. 

Next we add n extra Hasenbusch pseudofermions. The matrix Wn (for the heaviest 
pseudofermion) is identical to what we computed above; for the lighter pseudofermions, 

Wj=ML{mj)Ml{mj+i)-^M\^{m^+^){M\^{mj)MR{mj))-^MR{mj+i)M^^ 

(4.8) 
The theta function becomes 9{p'^ — "^^ \S,ij\'^{Xij — 1)). An analytic formula is unilluminat- 
ing. However, one discovers the following result: One of the eigenvalues of Wj is always 
less than unity. The other one is greater than unity, but as the number of pseudofermions 
increases, this eigenvalue falls to a value which is only slightly greater than unity. Thus 
the constraint on the lower end of the momentum integral relaxes and the tunneling rate 
rises toward unity, the deterministic result. (For example, for n = 8 and m = 0.05 the 
lower mass pseudofermions' largest eigenvalue ranges from 1.03 to 1.11. The highest pseud- 
ofermion has a largest eigenvalue of 2.5. For one pseudofermion the one relevant eigenvalue 
is equal to 3.8.) More pseudofermions enhance the tunneling rate. 

A graph of this behavior from a series of simulations is shown in Fig. y for m = 0.05. 
Most of the change happens with the first few pseudofermions. We are not sure whether 
this result has much practical use in QCD, since our algorithm is costly enough that many 
pseudofermions are simply too expensive. 
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Figure 6: Crossing rate from left to right at the discontinuity from a simulation of the model for 
TO = 0.05. 



The average (stochastic) AS is related to its variance as above. The eigenvalues of W 
also give us the value of AS at the crossing point: for n pseudofermions the result is 



ASL^R = Y,^{^m~l) 



(4.9) 



and 



ASn^L = Y.E(Kl-^)- 



(4.10) 



Thus with more pseudofermions, as A„j falls toward 1, the means of both distributions 
shrink. This is clearly illustrated with a plot of AS for eight pseudofermions, Fig. ^. 

The conclusion of this model study is that the asymmetric distributions we observe do 
not affect detailed balance. The value of AS is an indirect measure of the expected tun- 
neling rate, through the eigenvalues of the ratio of pseudofermion matrices at the crossing 
point. More pseudofermions should enhance the refraction rate. 

The model does not directly address the question of whether the chiral algorithm 
should have a higher tunneling rate than the nonchiral algorithm. However, in simulations 
restricted to the opposite chirality sector the shift in the spectrum is smaller than in the 
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Figure 7: AS" at the crossing in the model, with m = 0.1 and 8 pseudofermions. 

same chirality sector. This amounts to a larger value of m in the model. With bigger 
m the tunneling rate even with a single pseudofermion is closer to its deterministic value. 
This plus the use of the exact weight for the zero mode suggests suggests that the chiral 
algorithm will evolve more efficiently than the usual nonchiral HMC. 

5. Conclusion 

We have discussed a method to simulate an arbitrary number of flavors of overlap fermions 
in hybrid Monte Carlo. Besides removing the constraint on flavor number from HMC, the 
algorithm has some practical features. It avoids the problems in the simulation associated 
with zero-modes. In sectors of non-vanishing topology, this facilitates and significantly 
stabilizes the inversion. 

As far as we know, this method is only applicable to simulations with the overlap 
action, since it needs the Ginsparg- Wilson relation to relate the spectrum of D to that of 
D^D, that D'D commutes with 75, and that changes of topology can be observed from 
zero crossings of the kernel action. 
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